function PerformanceTop23

addpath('basic_functions')
addpath('../Data')

Nfactors = 3;

[Return,Z,IsF,Nt,Chars,CharsNames,~,~,~] = package_data_all_greeks;
Return = Return.ret_daily;

XXX = Z;
[N,T,L] = size(XXX);
RET = NaN(T,L);
mean_ret = NaN(L,1); tstat_ret = NaN(L,1);
for i = 1:L
    for t = 2:T
        ind = ~isnan(Return(:,t)) & ~isnan(XXX(:,t-1,i));
        if sum(ind) > 0
            xx = XXX(~~ind,t-1,i);
            nS = xx<=prctile(unique(xx),10);
            nL = xx>=prctile(unique(xx),90);

            rr = Return(~~ind,t);
            RET(t,i) = mean(rr(nL),'omitnan') - mean(rr(nS),'omitnan');
        end
    end
    mean_ret(i) = mean(RET(:,i),'omitnan');
    tstat_ret(i) = sqrt(T-2)*mean(RET(:,i),'omitnan')./std(RET(:,i),'omitnan');
end

%% set Z0 as the base
Z0 = Z;

%% only consider Top 23 characterisrtics
f = [find(strcmp(CharsNames,'RV -- IV')) find(strcmp(CharsNames,'MarketCap')) find(strcmp(CharsNames,'Delta')) ...
     find(strcmp(CharsNames,'RV')) find(strcmp(CharsNames,'Max10')) ...
     find(strcmp(CharsNames,'Assets')) find(strcmp(CharsNames,'Stock price')) ...
     find(strcmp(CharsNames,'Moneyness')) find(strcmp(CharsNames,'Option price'))  find(strcmp(CharsNames,'Gamma')) ...
     find(strcmp(CharsNames,'Vega')) find(strcmp(CharsNames,'Turnover'))  find(strcmp(CharsNames,'Stock return')) ...     
     find(strcmp(CharsNames,'Rkurt')) find(strcmp(CharsNames,'RV -- MFvol'))  find(strcmp(CharsNames,'Debt')) ...     
     find(strcmp(CharsNames,'MFvol')) find(strcmp(CharsNames,'Leverage'))  find(strcmp(CharsNames,'IV ATM')) ...
     find(strcmp(CharsNames,'Stock return11')) find(strcmp(CharsNames,'IV slope'))  find(strcmp(CharsNames,'Profitability'))  find(strcmp(CharsNames,'Cash to asset'))];

Z = Z(:,:,[f 47]);
[N,T,L] = size(Z);

% construct W and X matrices that will be used for IPCA
W = NaN(L,L,T);
X = NaN(L,T);
for t = 1:T-1
    W(:,:,t) = squeeze(Z(IsF(:,t+1),t,:))'*squeeze(Z(IsF(:,t+1),t,:)) / Nt(t+1);
    X(:,t+1) = squeeze(Z(IsF(:,t+1),t,:))'*Return(IsF(:,t+1),t+1) / Nt(t+1);
    Return(~IsF(:,t+1),t+1) = NaN;
end
clear t

[G, F] = IPCA(X,W,Nt,Nfactors);
r2 = IPCA_R2(G,F,Return,Z,IsF,X,W);
disp(['Total R_i ' num2str(r2.TotalR2_i) ', PricingError_i ' num2str(r2.RelPrcErr_i)])

AlphaR = NaN(N,T,'single');
for t = 2:T
    Betas = squeeze(Z(:,t-1,:))*G;
    AlphaR(:,t) = Return(:,t) - Betas*F(:,t);
end

%% now compute regular alpha
Z = Z0;
[N,T,L] = size(Z);

% construct W and X matrices that will be used for IPCA
W = NaN(L,L,T);
X = NaN(L,T);
for t = 1:T-1
    W(:,:,t) = squeeze(Z(IsF(:,t+1),t,:))'*squeeze(Z(IsF(:,t+1),t,:)) / Nt(t+1);
    X(:,t+1) = squeeze(Z(IsF(:,t+1),t,:))'*Return(IsF(:,t+1),t+1) / Nt(t+1);
    Return(~IsF(:,t+1),t+1) = NaN;
end
clear t;

[G, F] = IPCA(X,W,Nt,Nfactors);
Alpha = NaN(N,T,'single');

for t = 2:T
    Betas = squeeze(Z(:,t-1,:))*G;
    Alpha(:,t) = Return(:,t) - Betas*F(:,t);
end

%% now compute ALPHAs for portfolios
% compute strategies returns
UB = 90;
LB = 10;
strat_alpha = NaN(L,3);
strat_alphaR = NaN(L,3);

% load bootsrap distribution of gammas
load 'tables_outputs/GammaVariancesTop23' GammaBoot;
GammaBootR = GammaBoot;
load 'tables_outputs/GammaVariances' GammaBoot;

XXX = Z0;
ALPHA = NaN(T,L);
ALPHAR= NaN(T,L);
ZZ = NaN(T,L,L);
for i = 1:L
    for t = 2:T
        ind = ~isnan(Return(:,t)) & ~isnan(XXX(:,t-1,i));
        if sum(ind) > 0
            aa = Alpha(~~ind,t);
            aar = AlphaR(~~ind,t);

            zz = squeeze(Z(~~ind,t-1,:));
            xx = XXX(~~ind,t-1,i);

            if mean_ret(i,1) > 0
                nS = xx<=prctile(unique(xx),LB);
                nL = xx>=prctile(unique(xx),UB);
            else
                nS = xx>=prctile(unique(xx),UB);
                nL = xx<=prctile(unique(xx),LB);
            end

            ALPHA(t,i) = mean(aa(nL),'omitnan') - mean(aa(nS),'omitnan');
            ALPHAR(t,i) = mean(aar(nL),'omitnan') - mean(aar(nS),'omitnan');
            ZZ(t,:,i) = mean(zz(nL,:),'omitnan') - mean(zz(nS,:),'omitnan');
        end
    end

    % calculate variance of alpha
    ind = isfinite(ALPHA(:,i));
    vv = var(ALPHA(~~ind,i));
    for ell = 1:L
        for k1 = 1:Nfactors
            zf = ZZ(~~ind,ell,i).*F(k1,~~ind)';
            vv = vv + var(zf)*(var(GammaBoot(:,k1,ell)) + mean(zf)^2);
            for k2 = 1:Nfactors
                if k2~=k1
                    z2ff = ZZ(~~ind,ell,i).^2.*F(k1,~~ind)'.*F(k2,~~ind)';
                    cc = cov(GammaBoot(:,k1,ell),GammaBoot(:,k2,ell));
                    vv = vv + cc(1,2) * mean(z2ff);
                end
            end
        end
    end
    strat_alpha(i,1) = 100*mean(ALPHA(~~ind,i),'omitnan');
    strat_alpha(i,2) = sqrt(sum(ind))*mean(ALPHA(~~ind,i),'omitnan')./std(ALPHA(~~ind,i),'omitnan');
    strat_alpha(i,3) = sqrt(sum(ind))*mean(ALPHA(~~ind,i),'omitnan')./sqrt(vv);

    ind = isfinite(ALPHAR(:,i));
    vvr = var(ALPHAR(~~ind,i));
    for ell = 1:size(GammaBootR,3)
        for k1 = 1:Nfactors
            zf = ZZ(~~ind,ell,i).*F(k1,~~ind)';
            vvr = vvr + var(zf)*(var(GammaBootR(:,k1,ell)) + mean(zf)^2);
            for k2 = 1:Nfactors
                if k2~=k1
                    z2ff = ZZ(~~ind,ell,i).^2.*F(k1,~~ind)'.*F(k2,~~ind)';
                    cc = cov(GammaBootR(:,k1,ell),GammaBootR(:,k2,ell));
                    vvr = vvr + cc(1,2) * mean(z2ff);
                end
            end
        end
    end
    strat_alphaR(i,1) = 100*mean(ALPHAR(~~ind,i),'omitnan');
    strat_alphaR(i,2) = sqrt(sum(ind))*mean(ALPHAR(~~ind,i),'omitnan')./std(ALPHAR(~~ind,i),'omitnan');
    strat_alphaR(i,3) = sqrt(sum(ind))*mean(ALPHAR(~~ind,i),'omitnan')./sqrt(vvr);
end

YY = [strat_alpha(:,[1 3]) strat_alphaR(:,[1 3])];
bh_t(1:2,1) = BH(tstat2pval(YY(:,2)).*(YY(:,2)>0) + (1-tstat2pval(YY(:,2))).*(YY(:,2)<=0),0.05);
bh_t(3:4,1) = BH(tstat2pval(YY(:,4)).*(YY(:,4)>0) + (1-tstat2pval(YY(:,4))).*(YY(:,4)<=0),0.05);

ind = ones(47,1); ind(f,1) = 0;
YY = [strat_alpha(~~ind,[1 3]) strat_alphaR(~~ind,[1 3])];
YY = YY(1:end-1,:);

XX.ret = 100*abs(mean_ret(~~ind,1));
XX.alpha = strat_alphaR(~~ind,[1 3]);
XX.bh = bh_t(3);
save tables_outputs/PerformanceTop23 XX;

return